! **************************************************************************************
!  Rob Crain (2009), Centre for Astrophysics & Supercomputing, Swinburne University
! **************************************************************************************
! Code to sift through the Indices output by SelectParticles at each redshift 
! and identify the set of unique IDs. These particles are then traced to IC 
! coordinates and velocities, and then 'reverse Zel'dovich approximated' to get 
! their coordinates in the limit a -> 0. This particle set can then be used elsewhere
! to generate a mask for zoomed ICs.
! **************************************************************************************

program main
  use hdf5_wrapper
  use logs
  use param_module
  use myfiles
  use memory_module
  use header_vars
  use particle_data
  implicit none
  integer*4         :: isnap,ifile,itype,ipart
  character(len=16) :: csnap
  logical           :: is_multi

  ! Perform startup tasks and read the input files
  call startup()

  write(logunit,*)' indir    = ',trim(adjustl(indir ))
  write(logunit,*)' outdir   = ',trim(adjustl(outdir))
  write(logunit,*)' icdir    = ',trim(adjustl(icdir))
  write(logunit,*)' icfile   = ',trim(adjustl(icfile))
  write(logunit,*)' istart   = ',istart
  write(logunit,*)' iend     = ',iend  
  write(logunit,*)' hdf5_ics = ',hdf5_ics  
  
  ! Find total number of IDs - start with 'OriginalSelection' and then 
  ! trawl through 'SelectedParticles'
  write(logunit,*)' -finding NumIDs in original selection...'
  call flush(logunit)
  call read_original_selection_header()  
  NumPart_Total_AllSnaps(0:5) = NumSelected(0:5)

!  ! SPOOF
!  NumSelected(0) = 0
!  NumSelected(2) = 0
!  NumSelected(3) = 0
!  NumSelected(4) = 0
!  NumSelected(5) = 0

  write(logunit,*)' -NumSelected '
  write(logunit,*)'  Type0 = ',NumSelected(0)
  write(logunit,*)'  Type1 = ',NumSelected(1)
  write(logunit,*)'  Type2 = ',NumSelected(2)
  write(logunit,*)'  Type3 = ',NumSelected(3)
  write(logunit,*)'  Type4 = ',NumSelected(4)
  write(logunit,*)'  Type5 = ',NumSelected(5)
  call flush(logunit)

  write(logunit,*)' -finding NumIDs encountered in earlier snapshots...'
  call flush(logunit)
  do isnap = istart,iend
     call read_select_file_header(isnap)
     do itype=0,5 
        NumPart_Total_AllSnaps(itype) = NumPart_Total_AllSnaps(itype) + NumGathered(itype) 
     enddo
  enddo

  ! SPOOF
!  NumPart_Total_AllSnaps = NumSelected

  write(logunit,*)' -NumGathered from earlier snaps '
  write(logunit,*)'  Type0 = ',NumPart_Total_AllSnaps(0) - NumSelected(0)
  write(logunit,*)'  Type1 = ',NumPart_Total_AllSnaps(1) - NumSelected(1)
  write(logunit,*)'  Type2 = ',NumPart_Total_AllSnaps(2) - NumSelected(2)
  write(logunit,*)'  Type3 = ',NumPart_Total_AllSnaps(3) - NumSelected(3)
  write(logunit,*)'  Type4 = ',NumPart_Total_AllSnaps(4) - NumSelected(4)
  write(logunit,*)'  Type5 = ',NumPart_Total_AllSnaps(5) - NumSelected(5)
  call flush(logunit)

  write(logunit,*)' -Total number of ParticleIDs to be recovered:'
  write(logunit,'(A11,I9.9)')'  -Type0 = ',NumPart_Total_AllSnaps(0)
  write(logunit,'(A11,I9.9)')'  -Type1 = ',NumPart_Total_AllSnaps(1)
  write(logunit,'(A11,I9.9)')'  -Type2 = ',NumPart_Total_AllSnaps(2)
  write(logunit,'(A11,I9.9)')'  -Type3 = ',NumPart_Total_AllSnaps(3)
  write(logunit,'(A11,I9.9)')'  -Type4 = ',NumPart_Total_AllSnaps(4)
  write(logunit,'(A11,I9.9)')'  -Type5 = ',NumPart_Total_AllSnaps(5)
  call flush(logunit)

  ! Allocate memory
  call allocate_all_ids()

  ! Read the IDs
  call read_all_ids()

  ! Sort the IDs and discard duplicate IDs
  call sort_and_discard()

  ! Find fraction of particles that are from final snapshot
  do ipart=0,NumPart_Unique(0)-1
     if(UniqueSnapIDs%Type0(ipart) == iend) NumPart_Unique_Final(0) = NumPart_Unique_Final(0) + 1
  enddo
  do ipart=0,NumPart_Unique(1)-1
     if(UniqueSnapIDs%Type1(ipart) == iend) NumPart_Unique_Final(1) = NumPart_Unique_Final(1) + 1
  enddo
  do ipart=0,NumPart_Unique(2)-1
     if(UniqueSnapIDs%Type2(ipart) == iend) NumPart_Unique_Final(2) = NumPart_Unique_Final(2) + 1
  enddo
  do ipart=0,NumPart_Unique(3)-1
     if(UniqueSnapIDs%Type3(ipart) == iend) NumPart_Unique_Final(3) = NumPart_Unique_Final(3) + 1
  enddo
  do ipart=0,NumPart_Unique(4)-1
     if(UniqueSnapIDs%Type4(ipart) == iend) NumPart_Unique_Final(4) = NumPart_Unique_Final(4) + 1
  enddo
  do ipart=0,NumPart_Unique(5)-1
     if(UniqueSnapIDs%Type5(ipart) == iend) NumPart_Unique_Final(5) = NumPart_Unique_Final(5) + 1
  enddo

  ! Report...
  write(logunit,*)' -found (unique ids in all snaps : unique ids in final snap) '
  write(logunit,'(A11,I9.9,A3,I9.9)')'  -Type0 = ',NumPart_Unique(0),' , ',NumPart_Unique_Final(0)
  write(logunit,'(A11,I9.9,A3,I9.9)')'  -Type1 = ',NumPart_Unique(1),' , ',NumPart_Unique_Final(1)
  write(logunit,'(A11,I9.9,A3,I9.9)')'  -Type2 = ',NumPart_Unique(2),' , ',NumPart_Unique_Final(2)
  write(logunit,'(A11,I9.9,A3,I9.9)')'  -Type3 = ',NumPart_Unique(3),' , ',NumPart_Unique_Final(3)
  write(logunit,'(A11,I9.9,A3,I9.9)')'  -Type4 = ',NumPart_Unique(4),' , ',NumPart_Unique_Final(4)
  write(logunit,'(A11,I9.9,A3,I9.9)')'  -Type5 = ',NumPart_Unique(5),' , ',NumPart_Unique_Final(5)  
  call flush(logunit)

  ! Read initial conditions to get IDs, coordinates and velocities
  if(hdf5_ics == 0)then 
     call read_binary_ics()
  else
     call read_hdf5_ics()
  endif

  ! Match the selected IDs to the initial conditions
  call match_to_ics()

  ! Impose a 'reverse Zel'dovich' approximation to shift particles
  ! into a->0 state, where they are effectively force-free.
  call reverse_zeldovich()

  ! Write data to disc
  call write_to_hdf5()

end program main

 
